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Tensor network algorithms provide a suitable route for tackling real-time dependent problems in 
lattice gauge theories, enabling the investigation of out-of-equilibrium dynamics. We analyze a U(l) 
lattice gauge theory in (1+1) dimensions in the presence of dynamical matter for different mass and 
electric field couplings, a theory akin to quantum-electrodynamics in one-dimension, which displays 
string-breaking: the confining string between charges can spontaneously break during quench ex¬ 
periments, giving rise to charge-anticharge pairs according to the Schwinger mechanism. We study 
the real-time spreading of excitations in the system by means of electric field and particle fluctua¬ 
tions: we determine a dynamical state diagram for string breaking and quantitatively evaluate the 
time-scales for mass production. We also show that the time evolution of the quantum correlations 
can be detected via bipartite von Neumann entropies, thus demonstrating that the Schwinger mech¬ 
anism is tightly linked to entanglement spreading. To present the variety of possible applications of 
this simulation platform, we show how one could follow the real-time scattering processes between 
mesons and the creation of entanglement during scattering processes. Finally, we test the quality 
of quantum simulations of these dynamics, quantifying the role of possible imperfections in cold 
atoms, trapped ions, and superconducting circuit systems. Our results demonstrate how entan¬ 
glement properties can be used to deepen our understanding of basic phenomena in the real-time 
dynamics of gauge theories such as string breaking and collisions. 

PACS numbers: 05.10.Cc, 03.67.Ac, 11.15.Ha 


I. INTRODUCTION 

The mechanism of quark confinement stands as a key 
concept in our understanding of the fundamental inter¬ 
actions in high energy physics m- As a quark and 
an anti-quark are pulled apart, the energy stored in the 
gluon string connecting them grows linearly with dis¬ 
tance, due to the confining nature of strong nuclear forces 
described by quantum-chromodynamics (QCD). In gauge 
theories hosting dynamical charges, there exists a critical 
length scale at which the confining string breaks, creating 
particle-antiparticle pairs which reduce the energy den¬ 
sity in the string and give rise to the hadrons at the string 
edges [5]. 

The static properties of string breaking have been 
widely explored using a variety of lattice methods, 
wherein the effective string potential separating static 
charges can be extracted by the Polyakov or Wilson loops 
[MS]. However, the real-time dynamics of gauge theo¬ 
ries are usually biased by a severe sign problem, and as 
such cannot be accessed using lattice Montecarlo simu¬ 
lations mm- In this paper, we apply tensor network 
(TN) methods in order to study the real-time dynamics 
of a lattice gauge theory (LGT) with dynamical charges 
and quantum gauge degrees of freedom in one dimen¬ 
sional systems. In particular, we investigate the real¬ 
time string-breaking dynamics in Abelian U(l) LGTs in 
(l+l)d, which share with QCD the basic feature of con¬ 


finement. 

In recent years, efficient numerical methods based 
on TNs have found widespread application to the real¬ 
time dynamics of strongly correlated low-dimensional 
systems m- They are nowadays routinely used to tackle 
a variety of condensed matter and atomic physics prob¬ 
lems, such as the evaluation of spectral functions of 
low-dimensional magnets and the quench or controlled 
dynamics of ultra-cold quantum gases in optical lat¬ 
tices [T3H261 . While TN methods have been extensively 
applied to spin and Hubbard-type models, only recently 
it has been shown how TNs can provide an ideal plat¬ 
form for the investigation of gauge theories, for example 
in the study of the static properties of the Schwinger 
model, the low-energy mass excitation spectrum, and 
the dynamics of deconfinement in 2D Z 2 LGTs [27H37]. 
In particular, the quantum link model (QLM) formula¬ 
tion of LGT [38HF4Q] - gauge theories whose link Hilbert 
space is finite-dimensional - has been used to develop ef¬ 
ficient general-purpose TN algorithms to describe static 
and real-time dynamical properties of Abelian and non- 
Abelian LGTs including generic forms of matter fields, 
and to present different possible quantum simulator im¬ 
plementations on different platforms: trapped ions, cold 
atoms in optical lattices and circuit quantum electrody¬ 
namics (QED) f4H - [55) (see also [56j 57] for recent reviews 
and references therein). 

TN methods are based on variational tensor structure 
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FIG. 1: Left panel: Hilbert space and gauge invariant states of the QLM. (i) In the quantum link formulation, the gauge fields 
defined on the links are described by spins (in our case, S=l). (ii) Staggered fermions represent matter and antimatter fields 
on a lattice bipartition: on the even (odd) bipartition, a full (empty) site represents a particle (antiparticle), (in) Hilbert space 
and gauge invariant states of the QLM. The Gauss law, Eq. constrains the number of possible states at each lattice site. 
Notice that the Gauss law depends on the lattice site due to the staggered fermions. Middle panel: cartoon states for the 
different stages of the string breaking dynamics (see text). Here the leftmost site is an odd one. Right panel: sample simulation 
for the electric field dynamics when quenching an initial string state (B region) connecting two charges, and surrounded by the 
vacuum (A regions) for mn — 0 = g. Primary string breaking takes place in four stages (C-F), until an anti-string is created 
in place of the original string. The latter also decays durin g th e secondary string breaking. The shaded areas represent the 
wave-fronts estimated from entanglement entropies (see Sec. IV), which are directly related to the electric field evolution. 


ansatze for the many-body wave function of the quantum 
system of interest: the tensor structure is chosen to best 
accommodate some general system properties, e.g., di¬ 
mensionality, boundary conditions and symmetries, while 
a controlled approximation is introduced in such a way 
that one can interpolate between a mean field and an ex¬ 
act representation of the system. Being a wave function 
based method, one has direct access to all relevant infor¬ 
mation of the system itself, including quantum correla¬ 
tions, i.e., entanglement. In one-dimensional systems, an 
efficient tensor structure is given by the Matrix Product 
State (MPS) ansatz [12, 1T4] . defined as, 

IV’MPs) = Yi A a\ A a 2 ’ 02 ■ ■ • TV 1 1 5 )’ (!) 

a 

where the tensor A contains the variational parame¬ 
ters needed to describe the system wave-function, = 

1.. ..,d characterize the local Hilbert space, and = 

1.. .., m account for quantum correlations or entangle¬ 
ment (Schmidt rank) between different bipartitions of the 
lattice. Indeed, setting m — 1 results in a mean field de¬ 
scription, while any m > 1 allows for the description of 
correlated many-body states. Given the tensor structure, 
the tensor dimensions and coefficients are then optimized 
to efficiently and accurately describe the system proper¬ 
ties by means of algorithms that scale polynomially in the 
system size and m. Usually, these algorithms exploit the 
system Hamiltonian tensor structure, naturally arising 
from the few-body and local nature of the interactions, to 
efficiently describe the system ground state or low-lying 
eigenstates, or to follow the real- or imaginary-time evo¬ 
lution of the system itself. Indeed, in the TN approach, 
real- and imaginary-time evolution have no fundamental 
differences at the computational level as there is no sign 


problem, and limitations arise - only in some scenarios 
depending on the specific dynamics of interest as wit¬ 
nessed by the fast increasing literature appearing based 
on this approach m - from the amount of quantum cor¬ 
relations present in the system wave function. 

Here, we show how TN algorithms allow for the study 
of the real-time dynamics of LGTs, focusing on the 
string breaking in a paradigmatic confining theory - the 
Schwinger model [59l IBT] in a quantum link formula¬ 
tion. We characterize the real-time dynamics of the pri¬ 
mary and secondary string breaking and show that string 
breaking is intimately related to entanglement produc¬ 
tion in the system. A qualitative picture for string break¬ 
ing in our models, together with a typical result for our 
time-dependent simulations on a system of L = 100 lat¬ 
tice sites, is illustrated in Fig. [T] Even more importantly, 
our simulations allow us to track the entanglement evolu¬ 
tion during string breaking: as we will show below, string 
breaking and the so-called Schwinger mechanism are in¬ 
timately connected to entanglement propagation, which 
we address by evaluating the so called von Neumann en¬ 
tanglement entropy. Finally, we show that TN methods 
can be used to study scattering processes between bound 
states of LGTs: we develop a scheme to engineer me¬ 
son collisions m, and we show how, very surprisingly, 
the scattering not only reflects into an enhanced rate of 
particle-antiparticle creation, but it does affect drasti¬ 
cally the entanglement properties of the system, which 
stays significantly correlated well beyond the scattering 
time-window. 

The paper is structured as follows: in Sec. [IT] we 
present the system Hamiltonian and recall the TN algo¬ 
rithm we are using throughout this work. In Sec. |III| we 
present the results on string breaking and mass produc- 
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tion dynamics, including a discussion on how this phe¬ 
nomenon can be observed using a quantum simulation 
platform. In Sec. |IV| we show how entanglement follows 
the string breaking dynamics, providing a quantitative 
picture which underlines how entanglement entropies are 
directly tied to string-breaking. Finally, we present our 
result on scattering in Sec. [V] and draw a summary of 
our results in Sec. IVI1 

II. MODEL AND METHODS 
A. Model Hamiltonian: QED in (l+l)d 

QED in (1+1 )d, also known as the Schwinger 
model [59] . represents an ideal testing-ground for the 
benchmarking and development of new computational 
methods. Despite its relative simplicity, this model cap¬ 
tures fundamental aspects of gauge theories such as, e.g., 
the presence of a chiral symmetry undergoing sponta¬ 
neous symmetry breaking [59H6T1 [63H70] . Even more im¬ 
portantly, this theory, like QCD, displays confinement: 
differently from (3+l)d QED, in (l+l)d electrons and 
positrons are confined, and interact via a long-range po¬ 
tential which increases linearly with distance. Due to 
the large energy cost associated with the electric flux be¬ 
tween charges at large inter-charge distances, the electric 
flux string is unstable to particle-antiparticle creation, as 
in QCD, and string breaking takes place [58]. While this 
phenomenon, directly connected to the Schwinger mech¬ 
anism of mass production out of a vacuum, has long been 
debated, and notable insights have been provided using 
a variety of approximate methods, a full quantum me¬ 
chanical understanding of the complex real-time dynam¬ 
ics taking place during string breaking is lacking due to 
the computationally complexity of the many-body prob¬ 
lem [Ml EH [73]. 

In the Hamiltonian formulation, its dynamics are de¬ 
fined by the following form: 

h = -t^2[^lui x+ i'ipx+i+'ipi +1 u x , x +iip x 

X 

( 2 ) 

X X 

where are fermionic creation/annihilation opera¬ 

tors describing Kogut-Susskind (staggered) fermions (see 
Fig. [l), U x , are the gauge fields residing on the 
(x,x + 1) link, and we denote the strength of fermion¬ 
hopping (the kinetic energy of electrons and positrons) 
with £, the staggered mass of the fermions with m, and 
the electric coupling strength with g , where E XiX+ i is the 
electric-field operator. The gauge generator is given by 

(—\\ x — l 

Gx = ^Px^Px H - E x ^ x +i E x _i^ x + ~ 5 (3) 

so that [Ff, G x \ = 0 and all physical states |4/) satisfy the 
Gauss law G x \4t) = 0. While in the Wilson formulation 


U x , x +1 are parallel transporters acting on an infinite di¬ 
mensional Hilbert space, we focus here on a formulation 
based on QLMs, where the gauge fields are represented 
by spin-1 operators, U x>x+ 1 = S + x+1 , E x>x+1 = S^ x+1 
and, as such, act on a finite-dimensional link Hilbert 
space [15] . In particular, the electric field operator allows 
three possible states for the electric flux, constraining the 
physical states per site as described in Fig. [l] A detailed 
discussion of the quantum link formulation can be found 
in Ref. |38H4Q] . while in Ref. {311 it was shown how such 
quantum link formulation reproduces the phase diagram 
and quantum criticality of the continuum theory. 

B. String breaking and classical cartoon states 

String breaking is the process of cutting and short¬ 
ening the electric flux string that connects a particle- 
antiparticle pair by creating a new charge-anticharge 
pair [45]. Within our framework, a string consists of 
two charges creating non-zero electric flux between them. 
The charges are represented by appropriate boundary 
conditions (static charges) or by excitations of the mass 
field at the site of the fermion (dynamical charges). This 
is realized by an effective jump of a fermion from the 
site of one charge to the site of the second charge sat¬ 
isfying the Gauss law. The string of electric flux then 
follows from Gauss’ law. The charges force the links into 
a non-zero flux state, according to the configuration of 
the charges in either one direction or the other. Before 
embarking in a full quantum mechanical investigation of 
string breaking, we now discuss the classical (t = 0) static 
picture which provides a simple yet informative illustra¬ 
tion of the different stages of the string breaking mecha¬ 
nism. A set of cartoons of the classical states is provided 
in Fig. [l] 

Vacuum. - In the vacuum (A), neither mass nor electric 
field excitations are present. The vacuum energy is thus 
E 0 = -fm. 

String.- In the string state (B), two mass excitations 
are present at the boundaries, and all electric fields con¬ 
necting the two are also in the | + 1) state. The resulting 
string energy then takes the form 

n ^ 

^string — Eq = — (L — 1) + 2 171. (4) 

Pairs. - In the pairs state (C) all the masses are excited 
forming charge-anticharge pairs with an energy Flairs = 

^T+mL. 

Mesons. - In a confined phase, particle-antiparticle pair 
production can favor the establishment of a vacuum state 
between two static charges, with a pair of mesons at the 
boundary of the string (see (D)). The resulting energy is: 

F/mesons Eq — g + 4:171. (5) 

At the static level, string breaking takes place at a critical 
distance F c , above which the meson state is energetically 
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favored over the string state ( E str [ ng (L c ) = F mesons ): 


L c 


4 m 

= — +3. 

r 


(6) 


Antipairs and Antistring.- The antipair-state (E) and 
the antistring (F) denote the pair-state and the string 
with the electric flux having the opposite sign. 

At a dynamical level, string breaking takes place as 
a consequence of the Schwinger mechanism: in (l+l)d, 
the vacuum between two charges of opposite sign is un¬ 
stable against particle-antiparticle creation |58j mm, 
eventually leading to the electric field in the system flip¬ 
ping sign and to mass production. In real-time, this pro¬ 
cess develops following intermediate consecutive steps, 
and is schematically illustrated in Fig. [lj at very short 
timescales, particle-antiparticle creation takes place in 
the middle of the string, creating the pair state depicted 
in (C). Subsequently, the electric field in the center of the 
string relaxes to 0, and the external charges get screened, 
effectively forming mesons (D). At this point the process 
reverses, first establishing a state with anti-pairs (E), 
which finally decays into a string of oppositely signed 
electric field with respect to the initial state, an anti¬ 
string (F). 

String-breaking is a direct consequence of confinement: 
while in a deconfined phase such as QED in (3+l)d, it is 
possible to separate opposite charges at large distances 
due to Coulomb’s law, in a confined phase the corre¬ 
sponding electric field string breaks due to the effective 
potential between charges increasing as a function of dis¬ 
tance. However, the exact real-time dynamics of string 
breaking are inaccessible to classical simulations based on 
Montecarlo sampling due to a severe sign-problem. While 
classical-statistical approaches can provide remarkable 
insights in some parameter regimes isainHza such as 
small masses, unbiased numerical simulations for arbi¬ 
trary parameter regimes have been lacking. In the follow¬ 
ing, we present a systematic study of the string breaking 
dynamics using MPS techniques. 


C. Tensor networks for lattice gauge theories 

Tensor network algorithms are one of the paradigms 
for simulating quantum many-body systems in low- 
dimensions, both in and out of equilibrium, via a repre¬ 
sentation of the quantum state with a variational ansatz 
for wave-functions and/or density matrices [12]. 

For one-dimensional pure states, on which we focus 
here, the starting point is to consider a class of states 
of the form given in Eq. ([l]) and depicted in Fig. [2] with 
some fixed auxiliary dimension m and physical dimension 
d : for example, for spin one-half systems one has d = 2, 
while the dimension m depends on the states studied 
and on the desired accuracy. For ground states of one¬ 
dimensional gapped Hamiltonians, m is in general inde¬ 
pendent of the system size TV, while for critical systems, 
due to area-law logarithmic violations, the auxiliary bond 
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FIG. 2: Tensor network representations of a many-body quan¬ 
tum system, a) Any quantum state can be described via a 
tensor T of exponential dimension d N , where d is the dimen¬ 
sion of the local Hilbert space and N the number of sites, b) 
In a MPS representation, the wave function is characterized 
by N local tensors A, each one of dimensions m 2 d , where 
m is the maximum Schmidt rank allowed between different 
bipartitions, c) Gauge-invariant tensor network: the gauge 
invariant state can be represented via a tensor network state 
where a MPO imposes the gauge invariance while a varia¬ 
tional MPS accommodates for the detailed description of the 
wave function. 


dimension scales asmtx clog TV, where c is the central 
charge of the system (see Ref. 75] for a review). The 
case of out-of-equilibrium dynamics as considered here is 
more challenging, and no general picture is known, thus 
the bond dimension m has to be adapted to each specific 
case and convergence has to be checked comparing the 
results at increasing bond dimension ESE!- 

Global symmetries of the system of interest, such as 
particle number or global magnetization and even more 
complex non-Abelian global symmetries, can be embed¬ 
ded in the wave function ansatz given in Eq. 0 in 
an elegant way by promoting each tensor A^ 13 3+1 to a 
symmetry-sector-preserving tensor to a tensor that pre¬ 
serves every symmetry sector; that is, every index of 
the tensor is dressed with the corresponding symmetry 
charge number [T8H8T] . This symmetric formulation of 
tensor networks allows one to describe wave functions 
exactly and more efficiently with the desired quantum 
numbers, addressing each symmetry sector separately. 
Recently, it has also been shown that local symmetries 
- namely gauge symmetries - can be embedded in such 
description by generalizing the wave function ansatz to 
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a gauge invariant one We briefly recall 

such a construction in the rest of this Section, and tailor 
it to the model we will study in the rest of the paper. 

The Hilbert space of a gauge invariant system is given 
by the direct sum of every sector with different val¬ 
ues of the ’’Gauss law”, G x \^q LM ) = g x \^ qlm) and 
T-L = ®g x 'H gx . Dae to the gauge symmetry, there is 
no physical (e.g., gauge invariant) operator, including 
the Hamiltonian, that connects any two different sectors. 
Hence, starting with a quantum state partially defined 
by the values of the Gauss law, or a set of local (gauge) 
constants of motion, the time-evolved wave-function will 
remain in this sector under the action of the Hamilto¬ 
nian. In the QLM formulation of LGTs, a gauge invari¬ 
ant tensor description is immediately obtained if we use 
a ’’rishon” m or Schwinger representation of the gauge 
degrees of freedom (independent of the nature of the lo¬ 
cal continuous symmetry - Abelian or non-Abelian - and 
the dimensionality of the lattice m)- 

In our case, the U (1) QLM in (l+l)-d with a spin -1 per 
link, the spin operator or electric field E XjX +1 = S ^ +1 
is represented by a pair of Schwinger bosons E XjX +1 = 
\ ( n R,x+ 1 ~ n L,x) with a total occupation of two bosons 
per link, i.e. tir :X +i + til,x = 2. The first step to build 
an efficient tensor network representation of such states 
is to identify a local Hilbert space spanned by the states 
| a) defined on the tensor product of the fermionic matter 
field on a lattice site and of the rishon states on its left 
and right, that is \a) = |nR )X ,n^ jX , 7 i£, jX ). This allows 
for the projection of the state into the gauge invariant 
subspace, restricting the local bases only to the ’’physi¬ 
cal” states: in our model, this process results in the five 
gauge invariant states (uniquely identified in terms of the 
Schwinger rishons for even and odd lattice sites respec¬ 
tively and depicted in Figjl] left) given by [45] : 


odd: | 1 ) = 11 , 0 , 2 ), | 2 ) = | 2 , 0 , 1 ), 

| 3 ) = 11 , 1 , 1 ), 14 ) = 12 , 1 , 0 ), | 5 ) = | 0 , 1 , 2 ), 

even: | 1 ) = | 1 , 1 , 0 ), | 2 ) = | 0 , 1 , 1 ), 

| 3 ) = 11 , 0 , 1 ), 14 ) = | 2 , 0 , 0 ), | 5 ) = | 0 , 0 , 2 ). 


Operator (MPO) as depicted in Fig. [2j such that: 


Y&QLm) = A af 2 ■ 

A Pn — : 
• • ^a N 

fill ^7i >72 

Oi\ ,q:^ CX.2 , 0:2 
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(8) 


\a 1 a 2 ... ajv) ? • 


For the model we study, the MPO has bond dimension 
three (7 = 1,2,3) and the non-zero elements of the ten¬ 
sors B for the odd sites can be expressed as: 
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The tensors for the even sites can be computed in a sim¬ 
ilar way. Given the gauge invariant tensor structure in¬ 
troduced above, one can reformulate the standard algo¬ 
rithms for tensor networks to compute the ground state 
of the model, or, as we will do in the next Sections, 
compute the real time evolution of some initial state, 
e.g., by means of a Suzuki-Trotter decomposition of the 
time evolution operator acting on a pair of neighboring 
sites 123 E| . The parameters used in our calculation, to¬ 
gether with a discussion of the errors involved in the ap¬ 
proximations, are presented in Appendix [A] Finally, we 
mention that one can further simplify the tensor struc¬ 
ture in Eq. § to increase the algorithmic efficiency and 
that this construction is completely transparent with re¬ 
spect to the Abelian or non-Abelian nature of the gauge 
symmetry, resulting in a drastic simplification of numeri¬ 
cal analysis of non-Abelian lattice gauge theories ED [3l|. 


III. STRING BREAKING 


These are the only states allowed locally by gauge in¬ 
variance (notice that the local Hilbert space is different 
on even and odd sites as a consequence of the staggered 
nature of the fermions). A gauge invariant many-body 
state clearly exists in the tensor product of such local 
basis states. However, not all tensor product combina¬ 
tions are compatible with the spin representation S : with 
our choices, only the states with n^ x+ i + til :X = 2 are 
allowed. This additional constraint can be satisfied by 
modifying the tensor structure of the many-body state 
ansatz, that is by introducing a projector which acts 
on nearest-neighbor lattice sites and enforces the correct 
number of rishons on the link ED 134]. 

Hence, the gauge structure of systems with a local con¬ 
tinuous symmetry can be encoded in a Matrix Product 


In this section we present our results of the lattice 
gauge TN numerical simulations for the real-time dy¬ 
namics of string breaking: we focus on the time evolution 
of the electric and matter fields, quantitatively studying 
the properties of string breaking and of the Schwinger 
mechanism. In the following section we analyze the time 
evolution of quantum correlations during this process - 
an analysis enabled by the TN approach - and show that 
the two figures of merit are intimately related. 

The setup for our simulations is a dynamical string sur¬ 
rounded by the vacuum. The total lattice length TV = 100 
is chosen such as boundary effects are negligible for the 
timescales we investigate. Some typical results of the 
electric field time evolution are shown in the left column 
of Fig. [ 3 ] (Column A) for different values of the fermion 
mass m and electric field coupling strength g (hereafter 
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FIG. 3: Real-time evolution of a string of electric flux of length L = 20 embedded in a larger lattice (of length N §= 100) in 
the vacuum state (the initial cartoon state is sketched on the left side of the figure with the notation of Fig[Tj) . The electric 
flux real-time evolution (Column A) is shown with the evolution of the mass excitations (Column B) and the evolution of the 
bipartite von Neumann entropy (Column C) for mn — 0,g = 0 (Line 1), m = 0.25, g = 1.25 (Line 2), and m = 3 ,g = 3.5 (Line 
3). The von Neumann entropy is calculated using a bipartition of the system defined via a cut between lattice sites x and x + 1. 


we set t =» h = 1 and times are given in units of h/t). In 
the top row, the string freely breaks as for g = m = 0 no 
energy-cost is needed for such process to occur and the 
system evolves according to the free hopping Hamilto¬ 
nian: the electric field starts to oscillate between string 
and anti-string displaying primary, secondary (marked 
by vertical lines) and subsequent string breakings. More¬ 
over, the string propagates into the vacuum creating two 
diverging electric field excitation wavefronts. In the mid¬ 
dle row, representing the same scenario for non-zero mass 
and electric field coupling strengths, string breaking still 
occurs; however, the string evolution is damped and sta¬ 
bilizes the mean electric flux around zero. Finally, for 
large mass m and g , the large mass suppresses vacuum 
fluctuations while the large electric coupling prevents the 
occurrence of string breaking and thus the string does not 
decay. However, some dynamics still occur as we will see 
in more details in subsection |IIIB[ as fermion-antifermion 


pairs are created using the energy of the electric field and 
then annihilated, thus restoring the electric field excita¬ 
tions. 

To perform a quantitative analysis of string breaking, 
we repeat the same simulation for different values of m 
and g and analyze the mean electric field of the central 
six lattice sites of the string as a function of time. The 
results are reported in Fig. [4] where one can clearly see 
that different scenarios might occur: either the string 
breaks when the mean electric field drops below zero, or 
the electric flux remains positive throughout the whole 
evolution and no string breaking occurs. In the limit of 
m = g = 0 the oscillations show the highest amplitude, 
which is reduced by changing at least one of the two pa¬ 
rameters; as previously noticed, for high values of either 
system parameter the string never breaks. In the regime 
between these two extreme cases, we observe a third type 
of behavior: the electric flux tends to zero, however no 
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FIG. 4: Mean electric field of the central six lattice 
sites as a function of time r for the electric coupling 
g = 0.00,0.25,0.50,0.75,1.00,1.05,1.10,1.15,1.20,1.25,1.50 
(orange to light blue) for m — 0 (top) and mn — 0.25 (bot¬ 
tom). Primary (secondary) string breaking occurs when the 
mean electric field crosses the zero-line from positive (nega¬ 
tive) values. 

anti-string is formed: the oscillation is strongly damped 
and the system remains in the broken string state (com¬ 
pare with the middle row in Fig. [3]). These findings are 
summarized in Fig. |5} For g,ra ;<1 we observe the full 
string breaking dynamics with at least the partial for¬ 
mation of a string with opposite electric field flux after 
reaching the broken string state (red area). For g, m > 1 
string breaking was not observed and the dynamics are 
dominated by the interplay of the state of maximum pair 
creation and the original string (green area). Finally, the 
white area in between represents the region of parameters 
where we observe the string breaking with over-damped 
oscillations. 


A. String wavefront spreading 

During the string breaking process, a wavefront of elec¬ 
tric flux spreads outwards, as can be clearly seen in Fig. 
[3] (Panel A1 and A2). In this Section, we quantita¬ 
tively characterize the wavefront spreading by analyzing 
its spreading velocity and the oscillation intensity. 

In Fig. [6j we show the wavefront propagation as a func¬ 
tion of time for different electric coupling g for the zero 
mass case. The lower inset illustrates how we calculated 
such propagation: we follow the electric field excitation 
on one side of string by means of tracking the difference 
A E between the gauge field at some position x and the 
next nearest neighbor site x + 2. Further, we define the 



FIG. 5: State diagram of string breaking. The red area shows 
the parameters where the string breaks and evolves into a neg¬ 
ative string (E m ean < -0.15.Fmax). The white area represents 
the parameters where the mean electric field approaches zero 
and stays around that value (|F^ mea n| < 0.15F? max ). And fi¬ 
nally, the green area represents the parameters without string 
breaking (F/ mean ^ 0.15F/ max ). 

time when this difference displays a maximum s arrival 
of the wave front m- As can be seen from the lower 
inset of Fig. [6j where different colors represent differ¬ 
ent coordinates x, a wave-like propagation can be clearly 
identified. 

Following this scheme, we plot the position of the wave- 
front as a function of time in the main panel of Fig. [6] 
The result shows an approximatively linear spreading af¬ 
ter an initial transient time of about r ~ 2, with a ve¬ 
locity almost independent of the values of the electric 
coupling for g < 1. Increasing g starting from g = 1, the 
velocity increases as well, until for g > 1.5, where the 
results are inconclusive as increasing g leads to smaller 
wavefront amplitudes and consequently the errors bars 
prevent an accurate analysis to be carried out. However, 
for small g , the spreading velocity can be extracted di¬ 
rectly from Fig. [6j By fitting the values for r > 2/t and 
m — 0 we obtain a value of xe = 1.96 ± 0.02. 

Finally, in the upper inset of Fig. [6j we repeated this 
analysis for different masses and g — 0. The results 
clearly show that, for sufficiently large m/t > 4, the wave 
front spread velocity has an inverse linear dependence on 
the mass. All these results are in agreement with a theo¬ 
retical estimate obtained assuming the ends of the string 
as sources of excitations: in a quasi-free or weak coupled 
model, the speed is related with the band-width of the 
kinetic term resulting on an excitation spreading velocity 
proportional to v t h = 2/m. 

B. Schwinger mechanism 

During string breaking, the Schwinger model dynamics 
exhibit particle-antiparticle pair production as a conse¬ 
quence of the energy released from the external electric 
field string. This phenomenon is usually referred to as the 
Schwinger mechanism, and has been studied extensively 


















FIG. 6: Spreading of the wavefront of the electric held for 
m = 0 and g increasing from g — 0 (blue) to g = 1.5 (orange). 
A linear ht for r > 2 results in a velocity of ve = 1.96 =b 0.02 
(dashed green line). Inset (top left): Time needed for the 
wavefront to spread one lattice site as a function of m for 
g — 0. Inset (bottom right): Wavefront as a function of 
time evolving from the last site of the original string (x = 
Xi, blue) to the lattice site x = Xi + 16 (orange) for g — 0. 
The wavefront is calculated using the electric held difference 
= E{x + 2)-E(x). 
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FIG. 7: Time needed for the maximal mass production caused 
by the Schwinger mechanism as a function of the mass for 
<7 = 0,0.5,1,1.5, 2, 2.5, 3,3.5 (orange to light blue). 


since its hrst presentation in 1951 m- In the following, 
we provide a systematic investigation of the Schwinger 
mechanism in the context of the U(l) QLM. 

In our analysis, in the initial state dehning the string, 
the only two mass excitations present are the two dynam¬ 
ical charges which create the string itself. However, dur¬ 
ing the dynamics, the energy of the string is transformed 
into mass excitations. When the maximum mass produc¬ 
tion is reached, the particles start to annihilate, either to 



FIG. 8: Large-mass behavior of the mass-production time- 
scale for g = 0. The black line fits the data points with 
Tmax — (1.54d=0.02)/m. Inset: Log-log plot of the same data. 


break the string or to restore it. The time needed to reach 
the maximum mass production r max depends on the mass 
and electric coupling as shown in Fig. [7] It clearly 
displays two different behaviors, depending whether the 
electric coupling g is greater or less than one. For small-g, 
r max is maximal for m = 0 and decreases monotonically 
with m. This occurs, recalling the results of Fig. [5j in the 
regime where string breaking is observable. Indeed, these 
are the cases comparable to the top row in Fig. [3] On 
the contrary, for larger values of g, we observe the max¬ 
imum to occur for m > 0. In particular, the maximum 
is obtained at a point where the energy of the electric 
field matches approximately the energy needed to fully 
convert the string into particle pairs, i.e., m = g 2 / 4. 
This corresponds to the dynamics as displayed in the 
bottom row in Fig. [3j In the regime of high masses we 
can use second-order perturbation theory to estimate the 
mass-dependence of the mass-production time-scale ana¬ 
lytically, as the model can be approximated as decoupled 
double-wells potentials, resulting in r max = ^ 

We checked this approximation comparing the analytical 
estimate with the numerical results in Fig. [8] The best 
fit results in r max = (1.54 ± 0.02)/m, in good agreement 
with the theoretical prediction. 


C. Observability of string breaking in synthetic 
platforms 

Recently, the implementation of Abelian quantum link 
models has been envisaged on different platforms, such 
as ultra cold atom gases in optical lattices @511161152], 
trapped ions @H] , and circuit QED architectures HziEii. 
In this context, the possibility of investigating the real¬ 
time dynamics using TN methods provides an invaluable 
tool to benchmark experiment against theory in (l+l)d, 
and to address the role of possible imperfections in quan- 
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turn simulators. 

Unavoidable imperfections which will be present in any 
implementation can be detrimental to the observation of 
both ground state physics properties and real-time dy¬ 
namics. The former case has recently been investigated 
in the context of adiabatic state preparation, where it 
was shown how gauge-variant perturbations weakly af¬ 
fect the fidelity of the loading process [33] . Here, we fo¬ 
cus instead on the effect of gauge-invariant imperfections 
on the string-breaking dynamics. Different from gauge- 
variant terms, the role of gauge-invariant imperfections 
cannot be systematically addressed in an experiment us¬ 
ing, e.g., post-selection over the experimental data. 

Following the implementation schemes in Ref. 051071 
HB] , one of the most common form of gauge-invariant 
imperfections are nearest-neighbor interactions between 
matter and gauge fields: 

h, = t;J2 n *(s z x - t > x + s z x , x+1 ). (io) 

X 

with n x = 'ip'l'ipx- This form of the imperfection is usu¬ 
ally generated as a resonant term in perturbation the¬ 
ory to next-to-leading order with respect to t. While 
this implies t £, for realistic implementations the dif¬ 
ference in magnitudes between these two terms cannot 
be made arbitrary large: this will require small absolute 
energy scales, thus making other sources of more detri¬ 
mental imperfections such as temperature (for the cold 
atom implementation) and disorder (in the circuit QED 
implementation) dominant. At a qualitative level, this 
interaction term can freeze the system into a configura¬ 
tion where the matter fields remain pinned. This is due 
to the effective attraction generated by the nearby elec¬ 
tric field configuration. 

To estimate the effects of these imperfections on a 
quantum simulation of the dynamics considered in this 
work, we repeat the numerical simulations including re¬ 
alistic imperfections expected in a first generation of ex¬ 
periments. In the top row of Fig. [9] we show the string 
breaking evolution of the electric field, in the presence of 
i7/, with the same system parameters as used for Fig. [3j 
with imperfections of the order of 10% £ = 0.1. The re¬ 
sults including the imperfections still clearly exhibit the 
physics observed in the imperfection-free results. In gen¬ 
eral, even the quantitative dynamics is very well captured 
up to long-timescales, as illustrated in the bottom row 
of Fig. [9] The only exception are intermediate g and 
m values, where significant discrepancies (up to 50%) 
are observed for intermediate timescales (middle panel 
of the last row). In all other regimes, we could observe 
deviations up to a maximum of 15% caused by typical 
imperfections £ = 0.1. 

IV. ENTANGLEMENT DYNAMICS 

One of the key aspects of MPS-based methods is that 
they give full access to the wave-function during the real- 
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FIG. 9: Time evolution of the electric field including imperfec¬ 
tions of the type Hi = £ J2 x n *(S*- ljX + S* ym+1 ) with £ = 0.1 
at the system parameters m = 0, g = 0 (left), m — 0.25, 
g m 1.25 (center) and m = 3, g — 3.5 (right). The bot¬ 
tom row shows the difference to the result obtained without 
imperfection. 

time dynamics. By considering the dynamics of quantum 
correlations embodied by entanglement, this enables us 
to tackle the string-breaking problem from a fully com¬ 
plementary viewpoint with respect to the electric field 
and mass generation studies undertaken in the previous 
sections. The main question we want to address in this 
section is whether and to what extent entanglement plays 
a role in the string-breaking dynamics. 

In the last decade it has been shown that entangle¬ 
ment play a fundamental role in many-body quantum 
processes, from quantum critical phenomena, to quan¬ 
tum information theory and other fundamental aspects 
of quantum physics. Moreover, several aspects of quan¬ 
tum field theories, such as the static properties of con¬ 
formal field theories, have also been extensively studied 
using entanglement measures [74, 75] . Additionally, en¬ 
tanglement was shown to play a crucial role in the limits 
of classical simulations of quantum systems calling for 
the need of the development of quantum simulators to 
overcome such limitations [12!]. 

A common way to quantify entanglement for pure 
quantum states is by using the so-called Renyi entan¬ 
glement entropies, and, in particular, the von Neumann 
entropy. Given a pure state with the density matrix 
p = |'0)('^|, the entanglement entropy is given by the 
von Neumann entropy of the reduced density matrix 
p(x) = T v L -xP E5]: 

S{x) = -Tr {p{x) log 2 p{x)} . (11) 

S(x) is a measure of the entanglement of a bipartition at 
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the lattice site x. The entanglement entropy takes values 
between S(x) = 0 for a separable state (product state) 
and S(x) = log 2 d, with d being the size of the Hilbert 
space, for a maximally entangled state. 


A. Von Neumann entropy after string breaking 

In Fig. [3] (Column C) we plot the time evolution of the 
entanglement entropy for the three different cases con¬ 
sidered before (Panels 1-3): as it can be clearly seen, the 
entanglement evolution resembles the mass and electric 
field excitations, indicating how the two phenomena are 
related. Firstly, the vacuum fluctuations for small g and 
m generate not only mass and electric field fluctuations, 
but also a high amount of correlations. Moreover, the 
electric field wavefront is mimicked by the entanglement 
behavior, once again showing that not only do excita¬ 
tions propagate from the string, but also correlations do 
as well. Secondly, the string breaking process is clearly 
a two-steps process: first a correlated pair is created in 
between odd-even sites, and later in between the even- 
odd sites (blue-yellow checkerboard pattern inside the 
string in Panels Cl and 2). Finally, for large g and m, 
when the string does not break, the correlations behavior 
drastically changes as well. There, within the string, the 
correlations are built periodically only between odd-even 
sites, while in the vacuum the correlations are drastically 
suppressed (Panel C3). 

A transparent picture on how entanglement is gener¬ 
ated during the quench dynamics can be gathered by 
monitoring the entanglement growth at different points 
in space. In Fig. [TO] we show the entanglement time evo¬ 
lution for a set of system bipartitions and parameters, 
e.g., cutting the system in the middle of the string or in 
the vacuum. 

The solid lines are results obtained for m = 0 and 
g = 0, while the dashed lines correspond to the case 
m = 3 and g = 3.5. The orange line represents the en¬ 
tanglement growth for a partition at lattice site x = 20, 
therefore in the center of the lattice region starting in the 
vacuum. In the zero-mass case, one can see that the en¬ 
tanglement entropy grows almost linearly for most part 
of the evolution - signaling the vacuum instability against 
resonant particle-antiparticle pair production. Towards 
the end of the evolution, the linear growth breaks down 
as boundary effects start to play a role. The remaining 
solid lines in Fig. [To| represent the behavior of a partition 
between an even-odd lattice site (violet) and between an 
odd-even lattice site (blue) and display a more complex 
behavior. These results were obtained in the center of the 
string while it breaks up: the counter-phase oscillations 
indicate the competition of two states, together with an 
overall growth of the entanglement entropy, though not 
as fast as in the vacuum. As we have seen, in this regime 
the string breaking is a result of consecutive hopping pro¬ 
cesses: fermions hop on the lattice creating mass excita¬ 
tions followed by annihilations, which result in the string 



FIG. 10: Bipartite von Neumann entropy S(x) at m = 0, 
g — 0 (solid), m — 0.25, g — 1.25 (dot-dashed) and m — 
3, g — 3.5 (dashed). Partition at the center of the initial 
string (x — 51—blue, x = 50—violet) and in the vacuum (x = 
20—orange). 


breaking with two remaining dynamical mesons. This 
dynamical process goes on and after two hopping pro¬ 
cesses the anti-string is created. Fig. [To| shows that these 
dynamics are well captured by the oscillations of the en¬ 
tanglement entropy. Each maximum represents one hop¬ 
ping process: at the first maximum of the blue line the 
fermions are about to hop for the first time, while at the 
first maximum of the violet line the second hopping pro¬ 
cess is at its peak resulting in the string broken state. At 
around r « 4, the blue line does not display a maximum: 
this ’depleted region’ signals the anti-string state, where 
the last hopping event creating the anti-string is again 
the first hopping to break the negative string (and thus, 
it takes twice as long for the entropy to reach the next 
maximum). 

Differently from before, in the massive scenario 
(dashed lines), we see that the entanglement entropy for 
the vacuum stays close to zero as the large mass and 
electric coupling strongly suppress the particle-pair cre¬ 
ation that triggered the strong growth of the entropy in 
the previous case. Moreover, in the middle of the string 
the entanglement entropy is drastically affected: the blue 
dashed line initially behaves as the solid line in the mass¬ 
less case, reflecting the same mass excitation by pair cre¬ 
ation. However, the violet dashed line always remains 
close to zero as further evolution into the string broken 
state is energetically forbidden: the state evolves back 
into the string and the correlations between the even-odd 
sites cannot be created. The system is then oscillating 
between two almost degenerate states: the initial string 
state and the state made out of pairs. This results in the 
oscillating behavior of the entanglement entropy between 
zero and one. Finally, the third case with m = 0.25 and 
g = 1.25 (dot-dashed lines) lies between the two previous 
limiting cases: here the string breaks, but does not evolve 
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FIG. 11: Estimation of the spreading velocity v from the 
entanglement entropy S. The wavefront is defined as when 
the entropy drops below a certain threshold value with respect 
to the vacuum. The resulting fit using v — vs~ a/[log 10 (AS') — 
So] 2 gives us an estimate for the spreading velocity of us = 
(2.0 d= 0.2). Inset: Log-log plot of the same data adjusted by 
vs and So. 


into an anti-string. In the vacuum, the entanglement evo¬ 
lution is very similar to the first case (m = 0 = g) as the 
entropy grows almost linearly after a transient regime. 
However, the slope is reduced by the nonzero mass. The 
entanglement in the center of the string initially evolves 
as for the massless case, but after the first two hopping 
processes for r > 2, the oscillation turns into vacuum¬ 
like growth. This is a strong indication for non-periodic 
string breaking: the dynamics, although being unitary, 
resemble a dissipative process where the electric field en¬ 
ergy irreversibly disperses into the vacuum. This irre¬ 
versible behavior directly resembles what we observe in 
the electric field dynamics, which does not display any 
clear periodic signature. 


B. Entanglement propagation and wavefront 

Even more remarkably, the real-space particle creation 
and the entanglement dynamics are quantitatively tied. 
We concentrate on the signatures of the wavefront of 
the string imprinted on the evolution of the entangle¬ 
ment entropy. We consider the case m = g = 0 as it is 
characterized by the most pronounced wavefront, where 
the string with its slow entanglement growth is embed¬ 
ded in the fast growing vacuum (see Fig. [3j panel Cl). 
To characterize the entanglement spreading due to the 
wavefront, we exploit the fact that the entanglement en¬ 
tropy in the vacuum is constant in space even though it 
evolves in time. Therefore, far enough from both sides 
of the string there is a plateau of constant entropy much 
higher than the entropy in the middle of the string. Thus, 
to define the wavefront of entanglement spreading due 



FIG. 12: Scattering of two dynamical mesons using the system 
parameters m = 0, g — 8. The plot illustrates the time 
evolution of the electric field E(x) as a function of the position 
x. Lower panel: number of charges Nb = ^2 xE b n x i n the 
system during the evolution (blue: B = {1,..., 32}), number 
of particles present in the center (purple: B = 16), number 
of charges on either side of the center (coinciding lines red: 
B = {1,..., 15} and orange: B = {17,..., 32}). 


to the string, one can look for the lattice site at which 
the entropy plateau starts to decrease. We identify this 
point computing the difference of entropy between near¬ 
est neighbor bipartitions: tracking when this quantity 
become bigger than a given threshold allows us to char¬ 
acterize the entanglement wavefront spreading. 

In Fig. [Tl]we show the estimated spreading velocity for 
different values of the threshold: the limit for the thresh¬ 
old value going to zero gives an estimate of the spreading 
velocity. A power law fit results in a spreading velocity of 
vs = 2.0 ± 0.2 in very good agreement with the analytic 
estimate of vt — 2 and the result from the electric field of 
ve = 1.96 zb 0.02 demonstrating the intimate connection 
between entanglement and electric field spreading. 


V. MESON SCATTERING AND 
ENTANGLEMENT GENERATION 

Bound states are a fundamental component of gauge 
theories, and understanding their complex internal struc¬ 
ture is one of the most ambitious goals of computational 
physics [83]. In experiments, such internal structure is 
usually explored by colliding heavy ions, so that the 
energy released during the process can be released via 
particle-antiparticle creation. This makes the ab initio 
numerical simulations of such scattering processes chal- 
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FIG. 13: Scattering of two dynamical mesons. Main panel: 
Entanglement entropy S(x) using a bipartition between sites 
x and x + 1 as a function to time. After the scattering, the en¬ 
tropy significantly increases in the system: this is a direct sig¬ 
nature of enhanced quantum correlations. Right panel: S(x) 
at different times (see color bar), showing a clear plateau af¬ 
ter the collision, which enlarges as a function of time. The 
empty circles show the current position of the maxima of the 
electric-field which follow approximately the mesons center of 
mass. The dashed line represents S(x) generated by a single 
meson, while the green bar highlights the difference A S' to 
the entropy of the colliding mesons (difference between full 
and dashed line at r — 120, x% — 17). 


lenging, as MC simulations suffer from a severe sign prob¬ 
lem when tackling real-time dynamics. 

Here, we show how TN simulations allow to investi¬ 
gate meson scattering for the (l+l)d QLM using TNs. 
First, we present a general procedure to implement scat¬ 
tering processes between composite particles, and discuss 
the electric field dynamics after the collision. Then, we 
present results for the entanglement dynamics during and 
after the collisions showing that the meson collision is ac¬ 
companied by the creation of entanglement between the 
mesons themselves. Indeed, as we will show, the entan¬ 
glement is bounded by the propagation wavefronts of the 
particles after collision, and is characterized by a constant 
plateau of the entanglement entropy within the region. 


A. Electric field patterns during meson collisions 

In order to produce the scattering process, we shall 
start with two composite particles. Each particle is 
charge/anti-charge pair, and divided only by one link, 
namely a meson, with opposite momentum. For the two- 
meson problem, there is a simple picture in the strong 
coupling limit: the massless theory is a free massive bo¬ 
son (meson) theory that is expected to become weakly 
interacting once a small mass term is included. Hence, in 
the strong coupling region, a possible two-meson bound 
state is loosely bound, while in the weak coupling region 


it is tightly bound. 

We start the numerical simulation with the state repre¬ 
sented in the cartoon (D) in Fig. [lj two mesons separated 
by a vacuum state of ten sites, which can be straightfor¬ 
wardly be written in a simple, separable matrix prod¬ 
uct state with t = 0. We provide momentum to the 
mesons by adiabatically moving them from the bound¬ 
aries toward the center of the system: this is done by 
introducing a deep box-shaped potential that decouples 
the mesons from the rest of the system, with the only 
dynamics allowed being the oscillate between its position 
and a neighboring site. The box-potential is removed at 
time Ti = 17.4 when the meson is exactly at half oscil¬ 
lation: from that point on the mesons evolve freely with 
an effective momentum mostly in one direction, one to¬ 
wards the other and eventually colliding [91 ^. In order to 
avoid vacuum fluctuations during the process, we choose 
a large value of g = 8. Fig. [T2| shows an example of such 
a scattering process. In particular, it shows the absolute 
value of the electric field of two mesons approaching each 
other, colliding in the center and there parting again. 
While before the collision the meson are tightly bound, 
after the scattering process the electric field diffuses, and 
the corresponding wavefront has a significantly attenu¬ 
ated signal. In the lower panel of Fig. [l2j we monitor the 
time-evolution of the total particle number (blue), clearly 
indicating that this quantity is approximately conserved 
over the entire time-evolution, due to the large electric 
field strength, which suppresses particle-antiparticle cre¬ 
ation. 


B. Post-collision entanglement generation 

A classical-like picture of the scattering process pre¬ 
sented above, would have that two particles move against 
each other and then bounce back as there is not enough 
energy available to generate a more complex inelastic 
scattering. However, this picture is oversimplified, as this 
is a fully quantum process. Indeed one can, once more, 
monitor the quantum correlations generated during the 
scattering process. This is done in Fig. [l3j where we show 
the evolution of the bipartite entanglement entropy: one 
sees that entanglement is created and that it is mostly 
carried by the two mesons. In this parameter regime, the 
vacuum does not generate entanglement due to the very 
large value of g 2 . Studying the bipartite entanglement 
entropy for different bipartitions and times, one clearly 
sees that there are two regimes: before the scattering 
occurs, the entanglement is present only in the biparti¬ 
tion that cuts the mesons wave packets, indicating two 
electron-positron wave packets internally correlated but 
not sharing any quantum correlations among them. To 
the contrary, after the scattering, the two wave packets 
become highly correlated even when their two centers of 
mass are clearly separated (see Fig. 12 for times r > 100). 

The values of the entanglement entropy indicate that 
one ebit of quantum information has been created dur- 
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ing the scattering process. In the right panel of Fig. [l3j 
we present various cuts of the entanglement entropy pro¬ 
file taken at different times, together with a compari¬ 
son with the entanglement generated by a single meson 
moving through the lattice (dashed line). The difference 
of AS « 1 between the two cases (highlighted in Fig. 

, green bar) clearly shows that one additional ebit of 
entanglement has been generated during the scattering 
process: that is, that a singlet state has been created 
between the two indistinguishable mesons. The entropy 
has increased as the information on which process occurs 
(either mesons bouncing and moving back, or each one 
moving freely through the other one and keeping mov¬ 
ing with its own momentum) is completely lost. Notice 
that this is a direct consequence of the wave nature of 
the mesons: the classical case of two scattering particles 
in one dimension - with the constraint that no double 
occupancy might occur in a single matter site - would 
necessary results only in backscattering, since the two 
hard classical particles cannot pass through each other. 


VI. CONCLUSIONS 

We presented a detailed tensor network study on the 
real-time dynamics of a lattice gauge theory in the pres¬ 
ence of dynamical charges and quantum gauge fields. 
Within this approach, we have shown that one has di¬ 
rect access to all local quantities of interests - the time 
evolution of the mass, charge and gauge fields - and to the 
quantum correlation between bipartitions of the system 
by means of the von Neumann entropy. We investigated 
the primary and secondary string breaking in QED in 
(l+l)d represented by an S' = 1 quantum link model 
with staggered fermions. In this context, we studied the 
real-time evolution of the Schwinger mechanism leading 
to mass creation and annihilation by means of the in¬ 
terplay with the electric energy released by the string. 
We quantified key properties of these effects such as the 
mass production rate of the Schwinger mechanism and 
the velocity of the electric field spreading. Moreover, we 
unveiled the relation between string breaking dynamics 
and the entanglement spreading in the systems and we 
showed that it is possible to study scattering dynamics, 
characterizing not only mass and charge real-time evolu¬ 
tion but also the creation of quantum correlations among 
scattered particles. Finally, we showed that the pre¬ 
sented results can be in principle verified experimentally 
in possible future quantum simulations, as they appear 
to be robust with respect to the most common sources of 
gauge-invariant imperfections appearing in most of the 
proposed implementations. 

This work paves the way to systematic studies of real 
time dynamical phenomena in Abelian and non-Abelian 
LGTs in low dimensional systems. Indeed, the present 
approach can be straightforwardly generalized to more 
complex LGT and geometries, e.g., ladders or cylinders, 
and also can be studied in presence of an external en¬ 


vironment by means of, for example, the tensor network 
approach presented in [84]. Moreover, one can also study 
the continuum limit of the LGT (as already discussed 
for Wilson theories [30]): for QLMs, this can either be 
done using dimensional reduction [56], or by increasing 
the quantum link representation, similarly to what has 
been done in [29] . Notice that the gauge invariant ten¬ 
sor network formulation behaves favorably as the speed¬ 
up it grants scales as the number of rishons per link 
squared [34]. The unprecedented access to the entan¬ 
glement dynamics in LGTs will allow investigations on 
the role of quantum correlations in different contexts en¬ 
abling a deeper understanding of the quantum real-time 
dynamics of lattice gauge theories. Finally, exploiting 
the capability to prepare a wide class of complex states 
granted from recent developments in quantum optimal 
control of many-body quantum systems [85] |86], more 
complex dynamics could be investigated. One example 
would be to perform extensive studies of the scattering 
at different energies. 

Alongside, these methods can also be applied to study 
condensed matter systems to compute, e.g., response 
functions of antiferromagnets described by LGT (spin 
ices, Resonating Valence Bond models, etc.) which are 
very hard to evaluate using MC due to analytic continua¬ 
tion m- Moreover, real time dynamics of gauge theories 
is fundamentally interesting to ab initio investigations of 
scattering equilibration and pre-thermalization [55] . Fi¬ 
nally to benchmark and verify small quantum simula¬ 
tions whose proposal have recently appeared for different 
platforms (ranging from cold atoms to superconducting 
circuits and trapped ions) and that can be foreseen to be 
experimentally implemented in the next years [4TH55] . 

While writing this manuscript, we became aware of a 
related work on string breaking using TN methods [88] . 
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FIG. 14: Convergence tests of the electric field E for the cases 
m = 0, g = 0 (blue) and m — 0.25, g — 1.00 (orange). Left 
panel: Difference of the computed electric field from the offset 
E — E\[ m at t — S as a function of the inverse bond dimension 
X-The offset is obtained via a fit according toFoc(l/x) b + 
E \[ m . Inset: Truncation error during the time evolution for 
X = 200. Right panel: Difference of the computed electric 
field from the offset E — Eu m at r = 2 as a function of the 
time step St. The offset electric held is obtained via a fit 
according to E oc St b + Eu m . The exponent b is in perfect 
agreement with the expected value of 2. 


Appendix A: Convergence checks of the MPS 
simulations 


We investigated String-breaking in a one-dimensional 
system with up to 100 lattice sites. The simulation 
was performed using a matrix-product-state (MPS) al¬ 
gorithm using a second order Suzuki-Trotter decomposi¬ 
tion for the time-evolution. In this appendix we present 
the parameters we used for our simulations and provide 
a discussion on the relative errors. 

The main sources of numerical error in our calculations 
are the finite bond dimension of the MPS-state represen¬ 
tation and the Suzuki-Trotter time-step. As a bond di¬ 
mension we used a value up to x = 200, which ensures 
a truncation error on the corresponding wave function of 
maximum order 10“ 8 and 10“ 3 for r = 5 and 9 respec¬ 
tively for the string breaking calculations, and 10“ 5 for 
r = 250 for the scattering processes. Examples of the 
change of the truncation error during the evolution can 
be seen in the inset of the left panel in Fig. [l4j 

In order to make sure that these errors lead to small 
changes in our main observables, we tested the conver¬ 
gence of the mean electric field in the center of the chain 
(cf. Fig. [4]) for different values of x and typical system 
parameters. The results can be seen in the left panel 
of Fig. [14] where we plot the mean electric field at the 
end of the evolution (r = 8) as a function of the in¬ 
verse bond dimension 1/x- The mean electric field was 
subtracted by a fitted offset Eu m to allow a better com¬ 
parison between the two sets of system parameters. As 
we see from the plot, even for rather small bond dimen¬ 
sions x < 100 the change to the largest bond dimension 
used (xmax = 220) is in the order of 0.1 — 0.01%. For 
the bond dimension used in our simulations (x = 200) 



FIG. 15: Screening of external static charges after a time 
t = 400 using m = 0.05, g = 1, and x — 80 (blue squares). 
Two numerical lattice sites are combined to form one phys¬ 
ical lattice site due to our staggered fermions formulation. 
The dynamical charge density decays exponentially with the 
distance from the charges, the black line represents the the¬ 
oretical model with q{x) — aexp(— bx) — aexp(5x) where the 
decay parameter is b = g/y/n and the scale parameter is 
a — exp (5) — 1. 


the difference to the extrapolated correct value is in the 
order of E — E\[ m ~ 10 -5 . 

As said before, we used a second order Suzuki-Trotter 
decomposition with a time-step of Sr = 0.01 to simulate 
the time-evolution. In the right panel of Fig[l4]we report 
the convergence test obtained repeating the same simu¬ 
lation with different time-steps. As expected, we find 
a clear E — Eu m ~ Sr 2 dependence: This clean power- 
law behavior allows us to give a very good estimate of 
the correct value of E\[ m . In summary, the error from 
the Suzuki-Trotter decomposition using a time-step of 
St = 0.01 is of the order of E — Fi im ~ 10 -5 , that is, of 
the same order as the error introduced by the truncated 
bond dimension. 


Appendix B: Screening of Charges 

From the earliest discussion of the Schwinger model, 
the screening of static charges due to vacuum polariza¬ 
tion has attracted a lot of interest pH 151] . Recent semi- 
classical calculations have reproduced this process includ¬ 
ing also finite masses where it could be shown that the 
main dynamics for m < g is mostly independent from 
the mass . Here, we show that it is possible to study 
such process using quantum link models implemented in 
a tensor network algorithm, already with very small spin 
representation, e.g. S = 1. As the screening of charges 
can be observed only in asymptotically large times, the 
large correlations building up during the process have 
carefully to be kept under control, otherwise they would 
undermine the efficiency of our approach. Thus, here 
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we limit our system size to 20 lattice sites and encode 
the two static charges which build up the initial string 
in the boundary conditions on both sides of the lattice. 
In such a way, we avoid to simulate the vacuum regions 
where quantum correlations build up very quickly (see 
Fig.§ and are less interesting for the phenomena un¬ 
der study. Thus, we simulate the time evolution of the 
string of electric field among the charges with a small, 
but finite mass: after a long time compared to other 
timescales of the system (r = 400) when the largest fluc¬ 
tuations in the lattice damped out, and report the final 
charge distribution q x (averaged over the last r = 10 
to remove remaining oscillatory effects) in Fig. 15 The 
net results is an exponential decay of the charge density 
with the distance from the charges, as discussed in Ref. 
[89] where the authors showed that the screening also 
at m 0 is equal to that experienced by the massless 
Schwinger model, as q x = aexp(— bx) with a decay rate 
b = g/y/n ~ 0.5642 [51] . The parameter a = exp(6) — 1 is 
derived from the normalization condition Qx — 1- 


It can be clearly seen that our findings are compatible 
with the theoretical results from Ref. [89] (black solid 
line in Fig. 15). 


We stress that the goal of the presented analysis is 
to show that it is possible to perform a study of the 
screening of charges using tensor networks, not to per¬ 
form a thorough discussion which will be presented else¬ 
where. Indeed, an extensive numerical analysis is needed 
to individuate or rule out possible deviations of the mas¬ 
sive from the massless case in a full quantum mechan¬ 
ical model, the influence of the spin representation in 
the quantum link formulation, possible finite size ef¬ 
fect, and also more challenging regimes of parameters, 
e.g. m ~ 1 as shown in [89]. However, the agreement 
between the presented results and the semiclassical ap¬ 
proach, indicates that both analysis are capable of de¬ 
scribing the physics of the system and might be alterna¬ 
tive approaches which can complement and benchmark 
each other. 
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